Nonlinear screening and gas-liquid separation in suspensions of charged colloids 
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We calculate phase diagrams of charged colloidal spheres (valency Z and radius a) in a 1:1 
electrolyte from multi-centered nonlinear Poisson-Boltzmann theory. Our theory takes into account 
charge-renormalization of the colloidal interactions and volume terms due to many-body effects. 
For valencies as small as Z = 1 and as large as 10 4 we find a gas-liquid spinodal instability in the 
colloid-salt phase diagram provided ZXb/ci > 24 ± 1, where Ab is the Bjerrum length. 



Can like-charged colloids, dispersed in an aqueous sol- 
vent with monovalent cations and anions, spontaneously 
demix into a colloid-dilute "gas" phase and a colloid-dense 
"liquid" or "crystal" phase? For an index-matched solvent 
at room temperature, the answer to this question is no 
on the basis standard linear screening theory as formu- 
lated in the 1940s by Derjaguin, Landau, Verwey, and 
Overbeek (DLVO) [l|, |2|], but yes on the basis of some 
intriguing experimental observations in quasi-deionized 
suspensions of charged colloids 0, 3, @ • 

The classic DLVO theory, which is a corner stone of 
colloid science, predicts that pairs of colloids of radius 
a and charge — Ze at separation r interact through a 
screened-Coulomb potential 
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where the suspending medium is characterized by the 
temperature T, the Bjerrum length Ab = e 2 / {ek^T) with 
the dielectric constant e, and the Debye screening length 
At -1 . Here /cb is the Boltzmann constant and e the proton 
charge. The purely repulsive character of v(r) does not 
give rise to any cohesive energy to stabilize a dense liquid 
or crystal phase in coexistence with a much more dilute 
gas phase. Such cohesion effects, however, are observed 
in the experiments of Refs. [3, 0, 0J EJ- They include, 
for instance, compressed crystal lattice spacings indica- 
tive of gas-crystal coexistence 0], gas- like voids in an 
otherwise homogeneous liquid-like suspension long- 
lived metastable crystallites suggesting internal cohesion 
[H], and a (disputed) gas-liquid meniscus [f|. Explaining 
any of these phenomena requires cohesive energy, and 
the focus of much theoretical work has therefore been on 
extending DLVO theory to find "like-charge attraction". 

Important ingredients beyond DLVO theory are ion- 
ion correlations, nonlinear screening, and many-body ef- 
fects. For (dilute) 1:1 electrolytes correlations are consid- 
ered to be of only minor importance and will be ignored 
here, whereas the other two are important here. Non- 
linear screening has been studied extensively in the con- 
text of (spherical) cell models, where the geometry allows 
for straightforward numerical solutions of the nonlinear 
Poisson-Boltzmann (PB) equation @, S An im- 

portant nonlinear effect is quasi-condensation of ions onto 



a highly charged colloidal surface when ZX B /a > 5 — 10. 
As a consequence, the colloidal charge that appears in the 
prefactor of Eq.JT]) is renormalized from its bare value 
Z to a state point dependent Z* < Z, i.e., the inter- 
actions are reduced but remain repulsive. Interestingly, 
free-energy calculations on the basis of the nonlinear cell 
model show no sign of a gas-liquid spinodal instabil- 
ity [3| either. On the basis of these results, together 
with, e.g., formal proofs [ll| that nonlinear PB theory 
yields purely repulsive pair interactions, it is tempting to 
conclude that gas-liquid coexistence is impossible within 
(non)linear screening theory. 

However, there are also studies that do show cohe- 
sion and gas-liquid spinodal instabilities. Examples in- 
clude the early work by Langmuir [T3|, primitive model 

simulations of asymmetric electrolytes |l.'l eli m, pb 

calculations showing triplet attractions on top of pair- 
wise repulsions [16] in agreement with experiments [l3|, 
and the extended Debye-Hiickel theory and the boot- 
strap PB theory of Refs. Qj|. Interestingly, these 
systems have the explicit colloidal many-body character 
in common, as opposed to the cell geometry discussed 
above. Unfortunately, it is extremely time consuming 
and practically impossible to extend simulations such as 
those of Refs. [3, E3, E3| to realistic colloidal parame- 
ters (say, Z = 10 3 -10 4 with finite salt concentrations), 
or to calculate or measure effective n-body potentials 
with n > 4. Also attempts to study the full nonlin- 
ear many-body system, e.g., within the schemes as pro- 
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posed in Refs. 20, [2l|, [22|, turn out to be computation- 
ally difficult in the colloidal parameter regime. In fact, 
the colloidal parameter regime has so far mainly been 
studied as an explicit many-body system within linear- 
screening theory, where the many-body character ap- 
pears as a nontrivial density-dependence of pair inter- 
actions and as volume terms that can drive a gas-liquid 



transition [23, l2J, 123, 120, 122], albeit mostly in regimes 



where charge renormalization should have been taken 
into account thereby stabilizing the suspension (25| . 

The key question addressed in this Letter is whether 
the intriguing and sometimes hotly-debated experiments 
like those of Refs.[3|,|4|,|5|,l6| can be explained by hard-core 
repulsions and Coulomb interactions only. This ques- 
tion is answered by calculating phase diagrams of the 



2 



primitive model in the colloidal limit (Z> 1 and point- 
like cations and anions). We retain the best aspects of 
two well-established theories by combining the nonlinear 
screening effects of cell theory with the explicit many- 
body character of linear screening theory. 

We consider iV spherical colloids of hard-core radius 
a and fixed charge —Ze at positions {Ri} in a solvent 
of dielectric constant e and volume V at temperature 
T. The system is in osmotic contact with a salt reser- 
voir of monovalent point-like cations and anions at pair 
density 2c s , which gives rise to (yet unknown) ion con- 
centrations c± in the suspension. We only consider pair- 
wise Coulomb and hard-core potentials between any pair 
of particles, such that contact distances are 2a, a, and 
for a colloid-colloid, colloid-ion, and ion-ion pair, re- 
spectively. Thermodynamic properties and the phase di- 
agram of this system can be determined once we explic- 
itly know the semi-grand potential F(N, V, T, c s ), which 
describes the colloids canonically and the ions grand- 
canonically, and which is defined by the partition func- 
tion exp[-/3F] = (l/N\) J v dRi ■ ■ -dR N exp[-/3i2"{R}]. 
Here (3 — \/{k^T) and £f({R}) is the effective col- 
loid Hamiltonian. It was shown in Ref.[23] that H is 
the sum of the bare colloid-colloid interactions (hard- 
core and unscreened Coulomb) and the grand poten- 
tial of the inhomogeneous fluid of cations and anions in 
the external field of the fixed colloidal charge distribu- 
tion q(r) = — ZJ2iLi ^(l r ~ Ri\ ~ a)/47ra 2 . Within mean- 
field theory, the effective Hamiltonian can be written as 
a sum of entropic and electrostatic-energy terms [23] 
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where Has denotes the hard-core interactions between 
the colloids, and where all r-integrals are over the vol- 
ume outside the colloidal cores: |r — Rj| > a for all 
i = 1,...,N. The yet unknown quantities are the equi- 
librium density profiles p±(r) = c s exp[=F0(r)[ of the 
cations and the anions, and the (dimensionless) elec- 
tric potential 4>(r). Note that 4>{t) = in the reser- 
voir. These unknowns follow from Poisson's equation 
V 2 0(r) = — 47rAe[p+(r) — P-(r) + q(r)], which can be 
recast as the multi-centered non-linear PB equation 

V 2 0(r) = k 2 sinh</()(r) r outside cores (3a) 



n, • V0(r) = ZXb/o 2 



r = Ri + arij, 



(3b) 



where k 1 = 1/^87tAbc s is the reservoir's Debye screen- 
ing length, and where n.; is the unit surface normal of 
colloid i = 1, . . . N, 

In order to evaluate H of Eq.|f2|), we approximately 
solve Eqs. f3aj) - I|3bp as follows. We imagine each colloid 
i = 1, . . . , N in the center of a virtual cell of yet unknown 



radius b > a, and assume that the potential inside each 
cell is spherically symmetric and given by 4> c (\r— R^|) for 
a < |r — R;| < b. Denoting the net (yet-unknown) charge 
of the cell by Q, the cell-potential </> c (r) for a < r < b is 
the solution of the radially symmetric PB equation 



r 2 dr 
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The boundary value problem ([4]) is easy to solve nu- 
merically on a radial grid for given k, b, Ab, Z and Q. 
Outside the cells we retain the multi-centered character 
of 4>{v), but we exploit that it varies only weakly from 
some spatial constant (j> (provided b is large enough), 
such that 4>(r) — cf> is the small parameter in a linearized 
treatment of Eq. (|3aj) . The linearized multi-centered PB 
equation can be solved analytically [23], and in terms 
of a (yet unknown) effective colloidal charge —Z*e the 
resulting potential outside the cells is given by </>(r) = 
4> — tanh0 + YliLi 0l(l r — R-i | ) , with the Yukawa "or- 
bitals" 



^i(r) = -Z*\, 



cxp(Ka) exp(— nr) 
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(5) 



Following Ref.[23], we identify <j) with the Donnan po- 
tential of the renormalized system, such that sinh0 = 
-Z*n/(2c s exp[-T)/{l - rj)]), with n = N/V the colloid 
density and rj — 4ira 3 n/3 the packing fraction. Hence 
also c± — c s exp[q=0] and R 2 — 47tAb(c+ + c_) are known 
explicitly. Therefore, (j>(r) is completely determined in- 
side and outside the cells once Q, b and Z* are specified. 

How to determine these unknown quantities? In this 
Letter we calculate Q and Z* , for a fixed 6 to be discussed 
below, by imposing continuity of the (average) potential 
and its gradient at the cell boundaries. Choosing the 
origin at Ri, this implies that 

N 

<f> c (b) = 4>- tanh^ + 0i(&) + (Y,M\bni - R*|) 
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■' c (b) = Mb) + (ni • Vj2M\bni - R|)). (6b) 



The angular brackets in the right hand side of Eq. ([6ajl 
indicate an average involving the colloid-colloid radial 
distribution function g(R), and this term can be written 
as n JdIig(R)(f)i(\br — R|), and likewise for Eq. (|6b|) . In 
principle one could think of setting up a scheme to cal- 
culate g(R) and the effective hamiltonian ff({R}) self- 
consistently, yet here we are satisfied with the crude ap- 
proximation that g{R) = for R < 2b and 1 for R > 2b. 
This simplification allows for straightforward analytic ex- 
pressions for the bracketed terms in Eqs. (|6aj) and l|6b|) . 
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which can then be easily solved numerically for the two 
unknowns Q and Z* at fixed (Z, As, a, b, rj). 

The remaining problem is to choose the cell radius b. 
We checked that the largest physically reasonable choice, 
b = a?7 _1 / 3 , for which the cell has the volume of the 
Wigner-Seitz cell, yields essentially charge-neutral cells, 
Q/Z <C 1, such that (i) Z* is identical to Alexander's 
renormalized charge [7J, |8j , and (ii) no gas- liquid insta- 
bility is expected [9J. One also easily checks that the 
smallest possible choice, b — a, ignores any nonlinear- 
ity and leads essentially to the volume-term theories of 
Refs.EF ' 



which do predict gas-liquid spinodals. 
Our present choice for b is in between these two extremes, 
and is such that b is large enough for the region outside 
the cells to be properly described by linear screening, yet 
small enough to avoid significant overlaps of the cells. 
Moreover, since b is an unphysical parameter, the phase 
diagrams should not depend on b. All this is achieved, 
except perhaps in some extreme parameter regimes, by 
choosing b such that \4>c{b) — <j>\ — 8 (or b = a if Z\s/a is 
small enough), with fixed 8 ~ 1. This criterion for b leads 
to values of b/a changing monotonically from b/a — > 1 
(linear screening) for high packing fractions (77 > 0.5) or 
high salt concentrations (na > 5) to b/a > 10 (nonlinear 
screening) for dilute systems (rj <C 10~ 3 ) at extremely 
low salt concentrations {kcl < 0.1). 

Our numerical results show that Z* as defined by our 
procedure has all the characteristics of the renormal- 
ized charge as defined by Alexander i.e. Z* ~ Z if 
ZXb/o, < 5, if kcl is sufficiently large, or if 77 is high, 
and Z* < Z otherwise. For instance, for Z — 1000 and 
a = 100 As (which are of the order of the experiments 
of Refs.jlH), we find for t] < 10~ 2 that Z* ~ 700,900 
for na = 1,5, respectively, and Z* increases (in a non- 
monotonic fashion) to Z for increasing rj. These effects, 
which are independent of 8 (and hence of b) provided 
0.75 < 8 < 1.25, are also found in Alexander's traditional 
cell model [i^] and in experiments [29| . 

With Q and Z* determined from the continuity at the 
cell boundary (|6a|) and |6b"l) . and b from the criterion as 
described above, the potential <^>(r), and hence the ionic 
density profiles p±(r), are known explicitly, both inside 
and outside the cells. The Hamiltonian ([2]) can thus be 
evaluated explicitly. After tedious but straightforward 
calculation, in which the logarithmic terms are expanded 
to quadratic order outside the cells, one obtains 

N 



H({R}) = $ + v (\ R i ~ Rrfh Z *Vl - T(Ro), R) (7) 

i<j 

with a so-called volume term $ that does not depend 
on the coordinates of the colloids, and a sum of pair- 
wise screened-Coulomb interactions ^ with an effective 
charge and an effective screening parameter. Here we in- 
troduced the function T(x) = (1 + x)(l — exp[— 2x])/2x, 
such that the factor y/l — T(ko) in the effective charge 
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Figure 1: Typical r\ — na phase diagrams (insets) and crit- 
ical line Z\b/ci ~ 24 ± 1 (main figure), as determined over 
many decades of the parameters by the present theory, for 
suspensions of charged colloids (charge — Ze, radius a, pack- 
ing fraction n) in osmotic contact with a 1:1 electrolyte of 
reservoir screening constant k and Bjerrum length As. The 
critical line, which actually consists of three superimposed 
lines for S = 0.75, 1, 1.25 (see text), separates the strong- 
coupling regime (lower right) with a spinodal gas-liquid in- 
stability (dashed line) and a large density gap (grey area) 
at fluid-solid coexistence for low ko from the low-coupling 
regime (upper left part) without spinodal instability and with 
only a narrow density gap at fluid-solid coexistence. Criti- 
calit y a s found in the primitive model simulations of Refs.[l^| 
and [lj| are indicated by x and +, respectively, the simulated 
state point of Ref.[32] without an instability is indicated by □, 
and the experimental systems in which evidence for large den- 
sity gaps at gas-liquid and gas- crystal coexistence was found 
are indicated by ▲ [30], ■ [31], and • [||. 



is of order 1/2 in all but very extreme parameter regimes. 
The volume term $ in Eq.Q is given by 

j3§ f b 

= 4nc s n J drr {<j) c (r) sinh 4> c {r) — 2 cosh</> c (r)} 

2c_|_c_ 



-Z<f> c {a) 
b 3 



a 



_n 

1-7JC+ 
3" 



(8) 



77 



Z\ 



a=± 



111 - 



One easily checks that in the limit b — > a, for which 
Z* — » Z and 4> c {a) — > (f> — tanh</>+</Ji(a)-t-^^ 2 </>i(|ani — 
Ri|), one recovers essentially the linear screening theory 
of Ref.[23], and that the first term of the second line 
of Eq.® contains the Debye-Hiickel like — n 3 / 2 cohesive 
energy (per unit volume) that drives the spinodal gas- 
liquid instability at low salinity 23|, 2(| 2tJ. 

With the Hamiltonian explicitly given by Eq.Q, we 
can calculate phase diagrams from F = $ + -Fhsy , with 
Fhsy the free energy of a hard-sphere Yukawa fluid, 
which we calculate as in Ref.[23] by exploiting the Gibbs- 
Bogoliubov inequality. For fixed o/Ab and Z, we calcu- 
late phase diagrams in the rj-na representations by im- 
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posing equal osmotic pressure and chemical potential in 
the coexisting phases. The insets of FigQ]show two typ- 
ical phase diagrams, calculated with 6 = 1, where the 
upper left corner (Z — 10, a/A B = 100) only shows crys- 
tallization with a narrow density gap, and the opposite 
corner (Z — 2000, a/Ae = 10) exhibits a spinodal insta- 
bility at low na and hence a large density gap between the 
coexisting phases. From many such phase diagrams we 
constructed a curve in the (Z,o/Ab) plane of FigJU below 
which a spinodal instability is present in the rj-Ka plane. 
For 6 = 0.75, 1, and 1.25, these curves superimpose on a 
single line that is well approximated by ZXb /a = 24 ± 1 
over many decades of (Z, a/As). The independence of 6 
confirms again that the phase diagrams are independent 
of the artificial cell radius b (provided b is chosen sensibly, 
i.e. 0{6) = 1). The critical line from the linear screen- 
ing theory, b = a, satisfies ZXb /a — 6 — 7 (not shown) 
for Z > 10, where the weaker critical coupling indicates 
that charge renormalization has a stabilizing effect on the 
suspension in agreement with Ref.[25l|. 

FigQ] also shows the parameters of three colloidal sys- 
tems for which phase-instabilities have been observed 
experimentally (filled symbols in the figure); they are 
reasonably close to, yet systematically above, our pre- 
dicted critical line, by a factor of 1.5-2. Our critical 
line also shows good agreement with (estimates of) criti- 
cal points in the salt-free primitive model simulations of 
Refs. [3, Hi for 10 < Z < 80 (x) and Z = 2, 3 (+), 



respectively, although a systematically increasing devia- 
tion up to a factor ~ 3 (for Z = 80) is obscured by the 
double logarithmic scale of Fig.l. This deviation could 
perhaps be explained by the fact that for close-to-critical 
values of ZXb /a, instabilities (detached from the freez- 
ing line) first appear at na = but only for extremely 
dilute systems (77 <C 10~ 3 , decreasing with increasing 
Z), while the lowest density considered in Ref.[13l| is 
as "high" as r\ = 0.00125; the instability reaches this 
state point at a coupling that increases with Z 
Note also that for Z = 1, our theory predicts the crit- 
ical point at a/Ae ~ 0.048, which is close to the well- 
known critical point of the symmetric 1:1 electrolyte at 
a/Ae = 0.057 [34|. It is perhaps surprising, yet comfort- 
ing, that the critical point as predicted by the present 
theory, which is "designed" to deal with Z» 1, agrees 
quantitatively with the primitive model simulations of in 
the low- valency regime Z < 10. 

To conclude, we have constructed a theory for colloidal 
suspensions that interpolates between the well known 
limits of linear DLVO-type and non-linear cell-type PB 
theories, thereby combining the multi-centered charac- 
ter and the volume terms of the former with the charge 
renormalization of the latter. For high enough charges 
and small enough radii (ZXb/o, > 24±1, i.e. well into the 
charge renormalization regime) , we find spinodal instabil- 
ities at low ionic strengths. The theory directly extends 
the undisputed gas-liquid instability for asymmetric low- 



valency electrolytes (Z < 10) to the colloidal regime 
(Z > 1000), although the existing experimental evidence 
in the colloidal regime is at somewhat weaker couplings 
than required according to our predictions here. We are 
currently extending the present theory to take charge reg- 
ulation into account, in order to investigate whether the 
required coupling can be shifted towards experimentally 
realized systems. Alternatively it is interesting to con- 
sider the possibility that the experimentally determined 
"effective" colloid charge should not be identified with Z 
(as we did here) but with a renormalized charge [35j], as 
this would shift the experimental points of FigQ] closer 
to or beyond the critical line predicted here. 
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